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ABSTRACT 

Two-dimensional systems of linear hyperbolic equations are studied with 
regard to their behavior under a solution strategy that in alternate time- 
steps solves exactly the component one-dimensional operators. The initial 
data is a step function across an oblique discontinuity. The manner in which 
this discontinuity breaks up under repeated applications of the split operator 
is analyzed, and it is shown that the split solution will fail to match the 
true solution in any case where the two operators do not share all their 
eigenvectors. The special case of the fluid flow equations is analyzed in 
more detail, and it is shown that arbitrary initial data gives rise to "pseudo 
acoustic waves" and a non-physical stationary wave. The implications of these 
findings for the design of high-resolution computing schemes are discussed. 
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1. INTRODUCTION 


This paper is concerned with two-dimensional systems of time-dependent 
partial differential equations. We write these in their matrix form 


u H" Au + Bu « 0, 
~ t ~x ~y 


( 1 . 1 ) 


and restrict our attention to the hyperbolic case, characterized by the fact 
that the nxn matrix (Acos<j> + Bsin<j>) possesses a full set of inde- 
pendent eigenvectors for all angles <j> • Physically this corresponds to the 
fact that for any direction <J> it is possible to construct n different 
families of plane wave solutions travelling in that direction. 

A common technique for solving (1.1) numerically is that of operator 
splitting . The operation of advancing (1.1) through a time interval At is 
replaced, in the simplest case, by the following two operations. First the 
data is advanced by At by means of the one-dimensional operator 


u + 2Au = 0, (1.2) 

~t ~x 


and then the new solution is advanced through 



by means of 


u + 2Bu = 0. (1.3) 
~t ~y 

Repetition of this pair of operators N times is assumed to be equiva- 
lent to advancing (1.1) by NAt. The attraction of the strategy is that it 
allows use to be made of very accurate and sophisticated one-dimensional nu- 
merical schemes. The drawback is that the splitting itself introduces errors, 


proportional in the above case to At (although different sequences of split 

2 

operators reduce this to (At) ) even into a smooth flow [1]. 

Crandall and Majda [2] have pointed out additional errors that arise when 
the problem is nonlinear, albeit scalar. Singularities of the solution may be 
shocks, expansions, or contacts, corresponding to converging, diverging, or 
parallel characteristics. The errors occur when a singularity of one type is 
misread as being of another type by either of the one-dimensional operators. 

The error discussed in this paper is also a misreading, arising even at 
the linear level if there is more than one unknown, A and B are then con- 
stant matrices, and the problem is found whenever the eigenvectors of A are 
different from the eigenvectors of B, In that case, we will show the conse- 
quences of operator splitting when the initial discontinuity consists of a 
single oblique discontinuity and each of the one-dimensional split problems is 
solved exactly. The results do not therefore relate to any particular numeri- 
cal scheme, but can be used to explain certain anomalies that are observed in 
practice, and to suggest the kind of scheme that can be employed in given 
problems, A qualitative version of the analysis was given by Colella [3], but 
the present treatment reveals additional detail. 

We will study two specific cases. One of these is the 2x2 system 


u 

v 



0 li 

1 0 1 



(1-4) 


and the other is the 


4x4 system 


p 


u pa 2 0 0 


P 


v 0 p a 2 0 


P 

u 

_L 

1/p u 00 


u 

1 

0 v 0 0 


U 

V 

T 

0 0 u 0 


V 

+ 

1/p 0 v 0 


V 

s 

t 

0 0 0 u 


s 

X 

0 0 0 v 


s 


(1.5) 
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The first of these is an abstract system that serves to introduce the 
methodology; the second is a linear version of the compressible flow (Euler) 
equations, using pressure, velocity and entropy as unknowns. 

In Section 2 and 3 we will briefly recall the time evolution of an 
oblique inital discontinuity in the data (Riemann problem). In Sections 4 and 
5 we calculate how such a discontinuity evolves under the split operators when 
the governing equations are respectively (l.A) and (1.5). Section 6 contains 
numerical examples, and Section 7 discusses the implications. 


2. THE EXACT SOLUTION 

Consider initial data in which states u T , u_ occupy semi-infinite 

~L ~R 

regions separated by a straight line inclined at an angle 4> (Fig. 1(a)). 
Introducing coordinates x", y' as shown (1.1) becomes 

u + (Acosd> + Bsind>)u , + (Asind) - Bcosd>)u , = 0. (2.1) 

~t ~x ~y 

Solutions not depending on y' having the form 

u = f(x' - Xt)r 


can be found if X is an eigenvalue of (Acos<|> + Bsir*)>), and if r is 
the corresponding eigenvector. By a limiting process we can justify taking 
f to be the Heaviside function. To solve the proposed problem, project the 
given jump onto the eigenvectors 
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i=n 

~R"~L = l < 2 * 2) 

i=l 

In terras of Heaviside functions, the data is 

i=n 

u = u L + l a i (<(,)H(x')r 1 (<()), (2.3) 

and the solution at time t is (see Fig. lb) 

i=n 

u = ^ a i (<j>)H(x' - X i t)r ± (4> ) . (2.4) 

Our object is to discover how well the split operators reproduce this exact 
solution. 


3* THE OPERATOR-SPLIT SOLUTION 

During those periods when the equation being solved is (1.2), any jumps 
in the solution are eigenvectors of A. If this equation is used to evolve 
the initial data, the result is 


i=n 


u = u_ + Y a . (0)H(x'cos<j> - 2A (0)t)r (0) ■ 

~ ^ ii 0 1 1 ~ i 


(3.1) 


Similarly, if we use (1.3), we obtain 


i=n 


u = u^ + l Tr)H(x'sin<(> - 2X i (y lOOr^i- tt) < 


(3.2) 
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Now suppose that the initial discontinuity evolves for At according 


to (1*2). It will generate n distinct wavefronts (Fig. 2a shows the case n 
- 2). Each of these generates new data to be evolved for a further ^ At 
according to (1.3). At the end of the second interval we will have n^ waves 
in the system. After time kAt we would have n^k) waves. 


4. THE ABSTRACT SYSTEM 

Here the problem 

is defined by 


u ■+■ Au + Bu =0 
~ t ~x ~y 


i“i - ij _?i. »- 1? ji 


The eigenvalues of Acos<f> + Bsin<|> satisfy 


, cos<t> - X sin* | 

sin* -(cos* + X) 


i . e . j 


r » 1 . 


Devoting by r*(<J>) the eigenvectors for x « ±1, we obtain 


r (*) * 


In particular, 



-sin j <|) 


COS y * 

- 


» - 



1 J 


1 


COS j * 


sin 


.'0% - jOf 

— / r\\ _ 

f 1 1 

r 

vu; - |j | 


1 0 1 * 


(4.1) 


(4.2) 


(4.2) 


(4.3) 


(4.4) 


(4.5) 


and 


- 6 - 


r (j tt ) = 


-1//2 

+ ,1 x 

1//2 

1//T 

r (j tt) = 

1//2 


(4.6) 


'**’ ~h ~~ 

If these eigenvectors are denoted by r , r , r , r , we have the 

~x ~y ~y * 


identities 


(4.7) 


+ l , + -V 

r = — (r - r ), 

~ X /I ~ y ~ y 


r = — (r + r + ) , 
~ x ~y ~y 


+ 1 , - , +. 

~y ^ ~x ~x 7 ’ 


1 , - +. 
r = — (r - r ) . 

~y /j x x 


For example, we see that at the end of an x-step, each r + wave will 

~x 

give rise to r+ and r^ waves of amplitudes respectively 1//7 and 
-1/7 times that of the source wave. To develop the notation further we in- 
troduce translation operators C and S, corresponding to translation of a 
wavefront in the x'-direction by Atcos<f> and Atsin<|>. We use notation 
(=>) to mean "gives rise to" and derive from (4.7) 


r + =¥ — (Sr + - S ), 
~ x ~ y ~ y 


r =*> — (S _1 r“ + Sr + ), 
~ x ~y ~y 


r+ =$ — (C _1 r~ + Cr + ) , 
~y yj ~x ~x 


(4.8) 


r =s> — (C -1 r“ - Cr + ), 
~y ~x ~x 7 


This produces a calculus for tracing the history of a wave through any 
number of split operations. Thus, an r+ wave, operated on first in the x- 
and then the y-direction, evolves according to 


r 


+ 

x 


=» — (Sr + - S 1 r~) 

/2 y ~ y 

i [S(C"V + Cr+) - S _ 1 (C _ 1 r‘ - Cr+)] 

“ (SC + S -1 C)A I (sc -1 - S~ l C~ l )r\ 


Similarly, under the same operations 

r" =» i (SC - S -1 C)r + + \ (SC -1 + S _1 C _1 )r" . 
x 2 ~x 2 ~x 


* 

(4.9) 


(4.10) 


A powerful notation to describe the configuration of waves present after 

T 

a certain number of steps is to introduce a vector a = (a^,a 2 ) each of 

whose components is a sum of terras k a gS a C^. Such a term describes a wave 

of amplitude k^g at a location x' = asin<(( + gcos$. We use a^ to 

describe the r + waves and ao to describe the r - waves. Evolution 

~x * ~x 

through a pair of split operations from kAt to (k+l)At can then be 
described as 


*k+l 


Ma. . 
~k 


(4.11) 


where the matrix M can be obtained from (4.9), (4.10) as 

^ (SC + S -1 C) ^ (SC - S _1 C) 

y (SC -1 - S _1 C _1 ) y (SC -1 + S _1 C -1 ) 


M = 


(4.12) 
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The equation 


a. 

~k 



(4.13) 


+ — T 

is then a formal solution to the intial-value problem, if a^ = (a ,a ) 
where a + , a are the projections of the initial jump onto the eigenvec- 
tors r* . However, computing high powers of M would be cumbersome, and 
it is more convenient to use the recurrence relationship 


a k = trOOa^ - det(M)a k _ 2< (4.14) 

This follows from applying to a k _ 2 the identity 

M 2 - tr(M)M + det(M)I - 0 


which is the form taken for n = 2 by the Cayley-Hamilton theorem. Applying 
(4.14) to (4.12) gives 

a k - j (S + S -1 ) (C + C _1 )a k _ 1 - a^. (4.15) 

The iteration need not start with an isolated discontinuity, but may take a n 

~0 

arbitrary, a. = Ma rt . 

~ 1 ~0 

A particularly simple case that yields more general insights is to take 
♦ = ir, so that S = C = T (say) so that 

a k - \ <T 2 + 2 + T- 2 )a k .j - a^. 


(4.16) 
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Let c be the coefficient of T k in a, . Then (4.16) translates to 
~n ~k 


k+1 k 1 , k . k v k+1 
~n ~n 2 n-2 n+2^ ~n * 


(4.17) 


or, rearranging 


k+1 


c 

~n 


2c k + c k 1 
~n ~n 




2c 


+ 



(4.18) 


This shows that c^ behaves like a leapfrog approximation to the wave equa- 
tion 


-Stt 


2c = c , , . 

~XX ~X X 


(4.19) 


The approximation has an effective Courant number of 1//T . In (4.19) we 
see waves propagating to left and right with unit velocity as in the exact 
solution, but in (4.18) we will expect the phase error and oscillations that 
typify leapfrog schemes. 

The precise nature of the phase errors can be found by subjecting (4.18) 
to the usual Fourier analysis via the ansatz 


c k = g k exp(i9x'/At) 
~n 


(4.20) 


where g is a complex amplification factor, and 0 is the Fourier frequen- 
cy. Substituting (4.20) into (4.18) yields 


o 

g^ - 2gcos(0 siniji )cos(8 cos<|> ) + 1=0 


(4.21) 


or 


2 2 1/2 

g = cos(0sin<J>)cos(0sin<J>) ± i[l - cos (0sin<j))cos (0cos<|>)] . 


(4.22) 
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Since |g| 
is given by 


1 there is no dissipation. The phase shift per iteration 


tan~l — ~ cos ^ (8 sln<b )cos 2 (8 cos<|) ) ] 1 ^ 2 
an cos(8 sin<|> )cos(0cos<|> ) 


(4.22) 


which equals 0 whenever $ = Ntt/ 2. For other values of 4> , (4.22) can 

be expanded in 0 to give a phase velocity of 


sin 2 2<|> .2 

“TT - 9 


+ 0 ( 0 ^). 


(4.23) 


Since the operator-splitting solution has emerged as a nondissipative 
second-order approximation to the exact solution (for plane waves) standard 
results for such schemes apply. In particular, an initial discontinuity will 
be dispersed over an area proportional to and the accompanying 
overshoots will reach a maximum amplitude of 27% of the initial amplitude [4]. 


5. THE LINEARIZED EULER EQUATIONS 

In this section the problem 


u + Au + Bu = 0 
~t ~x ~y 

is defined by 


u = 

p 

u 

V 

> 

II 

u 

1/p 

0 

pa 2 0 0 

u 0 0 

0 u 0 

* 

w 

II 

v 

0 

1/p 

0 pa 2 0 
v 0 0 

0 v 0 


s 


0 

0 0 u 


0 

0 0 v 


(5.1) 


(5.2) 


The eigenvalues of (Acos<J> + Bsin<(>) are the roots of 
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det 


U - X 

COS<() /p 

sin<|>/p 

0 


pacos<() pasin<(> 0 

U - X 0 0 

0 U - X 0 

0 0 U - X 


(5.3) 


where U = ucos<|> + vsitv|> , i.e., 


(U - X ) 2 [ (U - X) 2 - a 2 ] = 0. 


(5.4) 


Corresponding to the two eigenvalues X = U are eigenvectors 


r 

~s 


0 

0 

0 

1 


* 


describing an entropy wave, and 


r 

~v 


0 

-sin0 
cosQ * 
0 


(5.5) 


(5.6) 


called here a vorticlty wave , across which the tangential component of 
velocity changes. 

Corresponding to the eigenvalue X = U + a is an eigenvector 


~a 


pa 

COS0 

sine 

0 


(5.7) 


which describes an acoustic wave . If 0 ranges over [— ir,ir] in (5.7), 
the eigenvalue X = U - a is accounted for as well. 

With no real loss of generality we can assume that the linearization is 
about a state of rest (both the exact and operator-split solutions are invari- 
ant under steady translations). In that case the entropy and vorticity waves 
are stationary, and the acoustic waves can move with the sound speed a in 
any direction (<j>). 
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To find the effect of the split operators we need to take <j> = o, y tt, 
yielding 



0 


0 


pa 1 


pa 


0 


0 

+ 

1 


— 

-1 

r = 

"-SX 

0 

• ^vx " 

1 

9 ~ax 

0 


* ~ax = 

0 
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0 



0 

1 

0 

1 

0 



pa 


pa 


0 


-1 

+ 


0 

— 

0 

1! 

CO 

0 

* 

In 

< 

ii 

0 

, r 

"-ay 


1 

* ~ay = 

-1 

1 

1 

l 

0 



0 


0 


( 5 . 8 ) 


( 5 . 9 ) 


The following identities determine how each time split wave will fission into 
waves travelling in the other direction 


r = r 
~sx ~sy 


1 r + " i 

r * ■=• [r - r ] 
~vx 2 ~ay ~ay 


+ 1 r + ^ i 

r =- T [r + r ] - r 

~ax 2 ~ay ~ay ~vy 


r = i [r + + r ] + r 

"-ax 2 "-ay "-ay "-vy 


r - r 
~sy sx 


1 r - + i 

r = tt [r - r ] 
"-vy 2 "-ax "-ax 


+ 1 r — + . 

r = tt [r +r ]+r 
~ay 2 ~ax ~ax ~vx 


( 5 . 10 ) 


( 5 . 11 ) 



From these results we can determine the evolution matrix for this prob- 
lem* Let a be a four-vector whose components describe the configurations 

of waves r + , r , r , and r respectively. Then some algebra reveals that 
~x ~x ~v ~s ° 


where 



(5.12) 


J (S + S -1 + 2)C J (S + S -1 - 2)C J ( s - S -1 )C 0 

j (S + S -1 - 2)C _1 j (S + S -1 + 2)C _1 j (S - S -1 )C _1 0 

j (s - s -1 ) T (s " s_1) 1 (s + s_1) 0 

0 0 0 1 

(5.13) 

and S, C represent displacements through aAtsin<(>, aAtcos<|>. 

The most obvious feature of (5.13) is that the entropy waves (at this 
linearized level) decouple completely from the rest of the waves, and are in 
fact represented exactly in the operator-split solution. It is easy to see 
that these statements are always true of any wave that is an eigenvector of 
both A and B (and hence of AcoS(|> + Bsin<J>). 

The matrix (5.13) has another property not observed in the earlier prob- 



lem (4.12). By treating S and C formally as real variables, we can 
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attempt to find eigenvalues of M. It is quite clear that there will be one 
eigenvalue A * 1, with eigenvector (0, 0, 0, 1) A corresponding to a con- 
figuration that is a steady solution of the system. The cubic that gives the 
other eigenvalues is 

(A 3 - I) - (A 2 - A)[y (S + S _1 ) + j (S + S -1 + 2)(C + C -1 )] - 0. (5.14) 


The existence of a second unit eigenvalue A = I implies another 
steady-state configuration. The corresponding eigenvector 


c l/2 (s l/2 . s ”l/2 ) 

-<f 1 / 2 (s 1/2 - s’ 1 ' 2 ) 

(C -1 / 2 - C 1 ^ 2 )(S 1 ^ 2 + S -1 / 2 ) 
0 


(5.15) 


itself exactly. (The other roots A 2 , A 3 of (5.14) are pseudo-operators, 
having no simple interpretations, although we may remark that A 2 A 3 - I. In 
this sense, they are inverse to each other and represent translations in 
opposing directions* We will describe the disturbances associated with them 
as ’’pseudo-acoustic” waves, and we shall see later that these waves obey a 
leapfrog-like recurrence relationship)* 

—11 2 

If, as an example, we take <j> * tan , so that C * S , the persis- 
tent configuration becomes 


s 3/2 _ s l/2 

S ' 3 ' 2 - s -‘ /2 

s -3/2 + 3-I/2 _ gl/ 2 . s -3/2 


(5.16) 


0 
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and this is sketched in Fig. 3a. All the waves that are radiated externally 
from this configuration will be self cancelling. The profiles of pressure and 
velocity are shown in Fig. 3b. Note that there is never any net change in any 
quantity across the persistent configuration (unless we add an entropy wave). 

For other values of <J> , the strengths of waves in the persistent config- 
uration remain unaltered, but they are found in different locations. The 
waves A, B, C, D are displaced from the origin by amounts 
■j aAt(cos<|> + sin<f>), aAt(cos<|> - sinj>) , -i- aAt(sin<f> - cos^i), - aAt(cos<|> + 
sin<f>), and occupy a region of total breadth aAt(|cos<j>| + sin|<Ji|). As <j> 
approaches Nir/2, however, the waves approach and cancel by pairs. 

Note that the persistent configuration could be multiplied by any poly- 
nomial in the shift operators, to yield another persistent configuration. 

We turn now to computing the evolution of arbitrary initial data. By the 
Cayley-Hamilton theorem we have 

M 4 - (D + 1)M 3 + 2DM 2 - (D + 1)M + I - 0 (5.17) 

where 

D = y (C + c" 1 ) + Y (S + S -1 ) + -J- (C + C _1 )(S + S -1 ) (5.18) 

so that 

a, = (D + l)a. - 2Da, + (D + l)a. , - a, . (5.19) 

~k ~k-l ~k-2 ~k-3 ~k-4 

is a valid, but tedious recurrence relationship. The trick is to separate a 
into its steady and unsteady parts, thus 


a = ap + Sr + n 

rsj At /wg /s) 
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I 


where 3 lies in £ , the complement of the space spanned by ^ and r g . 
The argument is formal, treating S and C as real numbers, but any inter- 
pretable statements reached in that way must be true. 

Introduce the eigenvector factorization of M 

M = R(A + A )L (5.20) 

P q 

- P + Q (5.21) 

where A^ = diag(l, 1 , 0 , 0 ) and A^ = diag( 0 , 0 , X^, X^) where X^, X^ 
do not need to be explicitly computed. We then have that 


9B = <£ s " p 3 " 0 


? 2 m 2 


Pr 


~s 


r 

~s 


Inserting (5.21) into (5.17) and using (5.22) gives 

(Q - I) 2 (Q 2 - (D - 1)Q + D 3 -O 


(5.22) 


and since Q - 1 is a non-singular, the part of the solution lying in £ 
obeys 

(Q 2 - (D - 1)Q + I ) 5 - 0 (5.23) 

or 


3 k - (° - D^k-i “ Sk-2 


(5.24) 
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Thus, the unsteady part of the operator-split Euler solution obeys a re- 
currence relationship very similar to that found for the abstract problem (cf. 
4.15). In the special case <j> = tt/4, we again find a leapf rog-type behavior 
for the discrete values 


3k - 2 3k-l + 3k- 2 - 4 ' 2 + T ' ! + T ~‘ + \ T " 2) 3k-r < 5 - 25 > 

This shows that the unsteady part obeys a numerical approximation to 
2 

= a Fourier analysis reveals a non-dissipative scheme; some 

lengthy algebra yields the phase speed in a general direction to be 

a(i + e 2 ) + o(e 4 ) 

which is only slightly in error. The principal disparity between the exact 
and operator split solutions is this. In both cases an arbitrary initial dis- 
continuity should project into a steady component (one that moves with the 
fluid) and an unsteady component. The chief failure of the split operators is 
that they do not correctly distinguish the two sorts of data (except for the 
entropy wave). 

In the exact solution, that part of the data that projects onto the 
acoustic eigenvectors is treated as unsteady. In the split solution, it is 
not easy to identify the unsteady part of arbitrary data, but it generally 
differs from the exact projection. We give here a partial analysis of the 
decomposition 

2 0 “ a 2 + 3 ~s + S 0 * 


(5.26) 
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This is facilitated by the observation that M can be symmetrized 


where 


M = E -1 X E 


(5.27) 


E = — 


1 _ 

/2 


c" 1/2 

c 1/2 

0 

0 


c 1/2 

c 1/2 

0 

0 

c _1/2 

-c _1/2 

0 

0 

E -1 = i_ 

c" 1/2 

-c _1/2 

0 

0 

0 

0 

1 

0 

/ 2 

0 

0 

2 

0 

0 

0 

0 

✓ 2 


0 

0 

0 

✓ 2 


(5.28) 


and 


X = 


1(C+ C 1 +2)(S+S 1 +2)-l i(C-C _1 )(S+S *+2) ^-(C 1/2 +C 1/2 )(S-S X ) 0 


!<C-cf 1 )(S+S -1 +2) ^(C+C _1 -2)(S+S *+2)+l i(C 1/2 -(f 1/2 )(S-S *) 0 


^(c 1/2 +c~ 1/2 )(s-s b ^(c 1/2 -c 1/2 )(s-s b y(S+S _1 ) 


(5.29) 


Since the eigenvectors of X will be orthogonal, a formal solution to the 
problem of determining a, 8 in (5.26) is 


a = 


»0 


r T E T E a 

o = ~s ~0 _ T 

3 TpTp ~s ~0 

r E E r 
~ s ~ s 


(5.30) 


(5.31) 
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and it is only the evaluation of (5.31) that presents any difficulty. Formal 
manipulations yield 


-[C- 1/2 (S 1/2 -S- 1/2 ),-C 1/2 ( S 1/2 -S- 1/2 )4(C- 1/2 -C 1/2 )(S 1/2 + S- 1/2 ),0] ; 

1/2 - 1/2 2 1/2 - 1/2 2 

8[1 -(£— S •) (5— £ ) ] 


(5.30) 


The denominator of this expression can be expanded as 


i T ,c 1 ' 2 + c- 1 ' 2 1 2 ",s 1 ' 2 + s- 1 ' 2 1 2 " 
a L \ -y 11 , ) 


(5.31) 


For (5.31) there is a coefficient of C^S^, say, given by a divergent 
infinite series. However, if (5.31) is multiplied by the numerator of (5.30), 
all coefficients are given by convergent, hence interpretable, series. For 
example, if the initial jump consists of an isolated r x wave, (ap = (1, 0, 
0, 0)) but inclined at <f> ■ 45° we have S ■ C * T and 


T" 1 - 1 ,T 1/2 + T -2 ^ 2 N 4n 

“ 5 — l ( 5 > • 

n=0 


(5.32) 


The largest coefficients are those of T“* and T®. We find 


c o = -°-i 


1 T (4n) ! 1 

' 8 n i 0 (2n) ! (2n+l) ! 16 n 

1 n r ^VlVlV^n 

~ "a L 12 2 3 


(5.33) 


(5.34) 


where (x) = T(x + n)/r(n) is the Pochammer symbol. Cancellation of terms 

n 
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leads to 

_ O 1 _ r 1 3 3 _ , 

C 0 = _C -1 " S 2 F 1 [ 7 » T’ 1 ; 11 

j r(|)r(i) 

8 r(f)r(£) 

_ i »>4> 

4 r4»4) 

This may be evaluated from the Reflection Formula 


r(z)r(l - z) = ir/sin(irz) 
as 

C Q = -C = j sin(ir/4) . (5.35) 


The numerical value of these coefficients is 0.177 which could be regarded as 
the amplitude of the persistent wave configuration generated by the particular 
initial data. The same amplitude would result from an initial jump propor- 
tional to r and inclined at <J> = 45°. 

If the initial jump is proportional to r yx and inclined at 45°, then 
(5.30) yields 


a 


T-T 1 / .T 1 ^ 2 +T -1 ^ 2 ^4n 

“XT' i c 2 > 

n=0 


(5.36) 


In this expression the coefficient of 


T° vanishes, but the amplitude may be 


represented by 
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16 
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( T ) n ( T ) n^7 ) n ( T^n 
16 n =° ( T ) n ( T ) n ( ? ) n ( T ) n 


i n 

= ^ I 


(5.37) 


i F [1 A 1*— 2*1] 

16 3 2 l 4 ’ 4 ’ ’ 2 ’ ’ J 


The generalized hyper geometric function appearing here can by summed by 
Watson's Theorem [5, p. 54], giving 

,2 


C = -C = i- 
L 1 -1 16 


r(^)r(|) 


r4>r(|) 


r4> 2 


-.2 


r 4 )r 4 ) 


which can be evaluated using the Reflection Formula as 


= sin 7T / 8 = 0.147. 


Results for arbitrary initial jumps can be found by superposing the 
special cases, but systematic conclusions would be hard to reach. The 

important point is that an appreciable fraction of the initial amplitude may 
be projected onto the persistent wave. This conclusion is reinforced by the 
numerical experiments in the following section. 
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6. NUMERICAL EXAMPLES 

Fig. 4 shows the distribution of pressure and vertical velocity though 
the wave system that results from 10 or 30 split time-steps being applied to 
data representing a shock inclined at <f> = tt/4, i.e., a^ = (y (1 + l//?), 
y (1 - 1//2), 1//2, 0) . Note that a persistent wave establishes itself very 
quickly at the center of the domain, where the initial discontinuity was 
placed. A typical "leapfrog" wave moves to the right in a manner that imi- 
tates rather loosely the exact solution (also shown). Some very weak disturb- 
ances move to the left. 

T" 

In Fig. 5 the initial data is a shear wave; a^ = (0, 1/2 , -1//2 , 0) . 
Here we see two "pseudo acoustic" waves moving outward, leaving a pronounced 
persistent pressure defect. The region between is filled with low amplitude 
noise, especially in the velocity. 

Fig. 6 demonstrates that the case <j> = x/4 is not untypical. Here we 
take <}> = tan *(y), computing first the data corresponding to a shock (a^ « 
(y (1 + 2r) , | (1 - 2r) , r, 0)^) where r^ = 1/5 and then due to a shear 
QJq = (0, 2r, -r, 0) t ). 

7. RELEVANCE OF STUDY TO PRACTICAL COMPUTATIONS 

We have seen that for discontinuous data, the solution of the split 
problem converges to the solution of the exact problem only in some weak 
sense, there being 0(1) errors even as At 0. That many practical 
calculations do successfully employ operator splitting is due to just one 
fact. Two successive wavefronts in the split solution are separated by at 

most aAt min[|sin$|, |cos<j>|] _< aAt. For explicit marching schemes this is 
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in turn less than the spatial mesh sizes Ax, Ay. In the frequency domain, 
the difference between the split solution and the true solution lies mostly 
above the folding frequency and so cannot be resolved by the majority of 
numerical schemes. 

This reasoning does not of course apply to implicit schemes, so any 
implicit method would have to contain "sufficient" dissipation. Nor does it 
apply to the explicit "large time-step" schemes developed by Harten [6] and 
LeVeque [7], In fact these latter schemes would be particularly inappropriate 
as they are "high resolution" schemes that aim at preserving high frequencies. 

Another scheme that performs very badly in split-operator applications is 
the Random Choice Method [3]. Again this is because the one-dimensional 
operator is just too good. The states generated in solving the one- 
dimensional Riemann problems are not lost in any averaging process, but have a 
chance to appear in the solution. Although only some of them will appear, 
there is a striking similarity between the results of the last section and the 
two-dimensional random choice results in [3]. It is to be expected that other 
attempts to create "very high resolution" schemes, such as Harten' s ACM [8] or 
Roe's Ultrabee [9] may experience similar difficulties. 

More constructively, it is worth noting that any scheme based on decompo- 
sition of the data into different waves could use a very-high-resolution 
operator-split method on those eigenvectors common to both operators, with a 
slightly less ambitious method employed on the remaining data. In the context 
of gas dynamics, we might use moderate resoluton on the acoustic waves. 
Shockwaves would nevertheless appear sharp, because they are nonlinear, "self- 
steepening" phenomena, and expansions would be treated well where they became 
smooth. The difficult (because linear) entropy waves could be sharpened with- 
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out creating any splitting problems. Even more valuable, perhaps, would be 
the application to chemically reactive or multi-phase flows. These always in- 
troduce new equations (for the species concentrations) that add more linear 
wavefields having the same eigenvector in all directions [10]. The ability to 
use very high resolution schemes on these fields will help to reduce numerical 
dissipation. 

The disturbance that seems most difficult to treat is the shear wave. 
Because it is linear it is prone to numerical diffusion, but a very high reso- 
lution splitting scheme would introduce into the solution the sort of noise 
displayed in Figs. 5(a-d) and 6(c,d). It may be necessary to control very 
carefully the degree of resolution applied to these waves, if they are not 
aligned with any coordinate direction. 


-25- 


REFERENCES 


[1] G. Strang, "On the construction and comparison of difference schemes," 
SIAM J. Numer. Anal. , Vol. 5, pp. 506-517, 1968. 

[2] M. Crandall and A. Majda, "The method of fractional steps for conserva- 
tion laws," Numer. Math. , Vol. 34, pp. 285-314, 1980. 

[3] P. Colella, "Glimm's method for gas dynamics, "SIAM J. Sci. Statist. 
Comput . , Vol. 3, pp. 77-110, 1982. 

[4] R. C. Y. Chin and G. W. Hedstrom, "A dispersion analysis for difference 

schemes: tables of generalized Airy functions," Math. Comp. , Vol. 32, 

pp. 1163-1170, 1978. 

[5] L. J. Slater, "Generalized hypergeometric functions," Cambridge, 1966. 

[6] A. Harten, "On a large time-step high-resolution scheme," ICASE Report 
No. 82-34, 1982. 

[7] R. LeVeque, "A large time-step generalization of Godunov's method for 
systems of conservation laws," SIAM J. Numer. Anal., Vol. 22, pp. 1051- 
1073, 1985. 

[8] A. Harten, "The artificial compression method for computation of shocks 

and contact discontinuities: III Self-adjusting hybrid schemes," Math. 

Comp., Vol. 32, pp. 363-389, 1978. 


-26- 


[9] P. L. Roe and M. J. Baines, "Asymptotic behavior of some nonlinear 

schemes for linear advection," in Notes on Numerical Fluid Mechanics , M. 
Pandolfi and R. Piva (ed.), pp. 283-290, Vieweg, 1984. 

[10] H. C. Yee and J. L. Shinn, "Semi— implicit and fully implicit shock- 

capturing methods for hyperbolic conservation laws with stiff source 
terms," AIAA Paper 87-1116, July 1987. 


- 27 - 



Fig. 1. Breakup of an oblique discontinuity. 
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